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摘 要 : 在 天 文 图 像 处 理 领域 ， 图 像 减 法 技术 得 到 了 广泛 的 应 用 。 由 于 不 同 的 大 气 
条 件 、 不 同 的 曝光 时 间 和 不 同 的 滤 光 片 等 原因 ， 两 幅 图 像 无 法 直接 进行 图 像 相 减 。 
本 文 提出 一 种 新 的 图 像 减 法 算法 ， 其 本 质 是 从 统计 学 意义 上 利用 相关 性 消除 两 幅 图 
像 之 间 具 有 相似 流量 分 布 的 部 分 ， 并 保留 具有 不 相似 流量 分 布 的 部 分 。 该 算法 执行 
速度 非常 快 ， 是 数值 稳定 和 局 部 独立 的 。 基 于 本 文 的 算法 ， 使 用 Python 语言 作为 接 
口 ，C 语言 作为 底层 实现 开发 了 一 套图 像 减 法 代码 。 使 用 该 算法 以 及 其 他 三 种 类 似 
的 算法 来 处 理 天 文 图 像 。 试 验 表明 ， 该 算法 能 够 在 很 短 的 时 间 内 找 出 两 幅 图 像 之 间 
的 差异 ， 并 快速 定位 运动 目标 ， 同 时 具有 良好 的 鲁 棒 性 和 位 置 测量 稳定 性 。 
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图 像 减 法 在 探测 未 知 星 源 、 寻 找 超新星 和 研究 微 引力 透镜 现象 方面 具有 重要 意 
义 。Tomaney 与 Crotts 首 次 在 频率 域 求解 卷 积 核 做 图 像 减法 趾 。 他 们 应 用 图 像 减法 来 
监测 微 引 力 透 镜 现 象 ， 并 处 理 了 M31 星 系 的 数据 。 最 优 的 图 像 减 法 由 Kochanski 等 人 
提出 己 ， 这 种 方法 是 一 种 非 线性 最 小 二 乘 拟 合 过 程 ， 即 使 对 于 稀疏 的 视 场 ， 并 且 只 
拟 合 一 小 部 分 像素 时 ， 计 算 时 间 也 很 长 。Alard 与 Lupton 也 提出 了 一 种 最 优 图 像 减 法 
算法 中， 该 算法 在 图 像 空 间 域 直接 求解 卷 积 核 ， 其 实质 是 基于 标准 线性 最 小 二 乘法 
求解 卷 积 核 。 这 种 方法 将 卷 积 核 分 解 为 多 个 高 斯 函数 乘 以 多 项 式 的 形式 。Alard 使 用 
空间 变化 卷 积 核 并 优化 了 这 种 方法 的 计算 策略 外 。Wozniak 使 用 Alard 的 这 种 算法 处 
理 了 380 GB 的 OGLE-II bulge 微 引力 透镜 数据 外。 

在 以 上 方法 中 ，Alard 与 Lupton 提 出 的 方法 最 为 经 典 。 但 是 ，Alard 与 Lupton 方 法 
要 求 使 用 者 定义 使 用 的 高 斯 基 函 数 的 数量 、 它 们 相关 的 sigma 值 和 多 项 式 的 阶 数 ,， 这 
一 系列 参数 的 选择 无 颖 会 使 用 户 感到 困惑 ， 且 需要 进行 大 量 实验 才能 获得 特定 数据 
集 的 最 优 参数 负 。 后 来 ， 基 于 Alard 与 Lupton 的 思想 ， 人 们 开发 出 了 多 个 类 似 的 图 像 
减法 技术 。 例 如 ，Bramich 将 卷 积 核 视 为 离散 像素 阵列 ， 并 直接 使 用 线性 最 小 二 乘法 
求解 卷 积 核 的 像素 值 ， 卷 积 核 的 大 小 成 为 唯一 需要 调整 的 参数 时 。Yuan 与 Akerlof 使 
用 交叉 卷 积 来 避免 对 高 质量 参考 图 像 进行 比较 的 要 求人 中。 这 两 种 方法 都 称 为 去 卷 积 
方法 , 是 在 Alard 与 Lupton 方 法 的 基础 上 发 展 起 来 的 ,这 些 方法 的 缺点 是 需要 调 参 中 。 
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Bramich 方 法 的 计算 速度 非常 依赖 于 卷 积 核 的 大 小 ,， 随 着 卷 积 核 的 增 大 ,计算 时 间 变 
得 难以 忍受 。 此 外 ， 这 些 算法 不 是 局 部 独立 的 ， 处 理 整 个 大 图 像 会 存在 困难 ， 通 党 
需要 将 图 像 分 成 几 个 较 小 的 部 分 分 别 进行 处 理 钻 ， 而 且 它 们 在 数值 计算 上 是 不 稳定 
的 站 。 输 入 图 像 或 参考 图 像 上 的 坏 像素 或 饱和 像素 会 对 一 组 在 卷 积 核 区 域内 的 像素 
造成 污染 。 因 此 ， 这 些 方法 通常 在 做 图 像 减法 之 前 需要 剔除 坏 像 素 和 饱和 像素 ， 并 
且 应 避免 使 用 过 大 的 卷 积 核 申 ,在 某 种 程度 上 , 这 些 操 作 增加 了 图 像 减 法 的 复杂 性 ， 
并 可 能 导致 一 些 有 意义 的 像素 被 忽略 。 直 到 Zackay 等 人 提出 了 一 种 基于 基本 统计 原 
理 的 算法 馈 ， 该 算法 颠覆 了 传统 的 基于 去 卷 积 的 图 像 减法 算法 的 原理 ， 是 局 部 独立 
且 数 值 稳 定 的 。 

本 文 也 基于 基本 的 统计 原理 提出 一 种 较 少 参数 的 可 蔡 代 方法 ， 依 据 输 入 图 像 和 
参考 图 像 之 间 的 流量 分 布 情况 来 做 图 像 减法 。 即 利用 两 幅 图 像 同 一 位 置 的 两 个 小 区 
域 之 间 的 流量 分 布 的 相关 性 来 消除 两 幅 图 像 间 相似 的 部 分 并 保留 不 相似 的 部 分 。 基 
于 本 文 的 算法 ， 开 发 了 一 套 以 Python 语言 为 接口 ，C 语 言 为 底层 实现 的 图 像 减 法 代 
码 。 该 算法 能 在 极 短 的 时 间 内 搜索 出 运动 目标 ,从 而 克服 人 眼 搜 索 运 动 目标 的 困难 ， 
本 文 称 之 为 消除 相似 和 保留 差异 ， 简 称 esapd(Eliminate similarity and preserve 
difference). 

1 方法 

为 了 搜寻 运动 天 体 ， 本 文 提 出 了 一 种 基于 相关 性 的 图 像 减 法 技术 ， 其 核心 在 于 
使 用 皮尔 逊 相 关系 数 求 出 两 幅 图 像 的 两 个 相对 应 小 区 域 的 相关 性 强 弱 ， 当 相关 性 较 
强 时 ， 缩 小 两 幅 图 像 之 间 的 差异 ， 当 相关 性 较 弱 时 ， 保 留 两 幅 图 像 之 间 的 差异 。 

如 图 1 AIAN, FA (c) 是 子 图 (a) AFA (b) 两 幅 图 像 的 差异 图 。 即 两 幅 图 
像 对 齐 ， 分 别 减 去 自己 的 背景 和 归 一 化 《〈 归 一 化 在 下 文 介绍 ) 后 直接 相 减 得 到 的 图 
像 。 子 图 (d) 是 子 图 (a) AFA (b) 两 幅 图 像 的 相关 图 (相关 图 的 获取 在 下 文 介 
绍 ) 。 实 验 发 现 ， 在 共同 背景 区 域 ， 相 关 图 的 系数 主要 分 布 在 -0.2 到 0.2 之 间 ( 相 
关 图 的 灰色 区 域 ); 在 共同 非 背 景区 域 , 相关 图 的 系数 主要 分 布 在 0.8 到 1.0 之 间 ( 相 
关 图 的 白色 区 域 。 运 动 目标 所 在 区 域 ( 非 重 车 情况 ， 即 运动 目标 没有 和 其 他 天 体 
HAUL. EARME PIC) 的 相关 性 也 是 较 弱 的 。 根 据 这 一 发 现 ， 可 以 缩 
小 相关 性 较 强 的 区 域 的 像素 灰 度 值 的 差异 ， 保 留 相关 性 较 弱 的 区 域 的 像素 灰 度 值 的 
差异 ， 即 可 突出 运动 目标 。 运 动 目标 的 位 置 如 子 图 (c) 的 红色 圆圈 位 置 所 示 ， 白 色 
目标 是 子 图 (a) 中 的 运动 目标 ， 黑 色目 标 是 子 图 b) 中 的 运动 目标 。 

为 获得 相关 图 ， 使 用 皮尔 偿 相 关系 数 ， 其 计算 方式 如 下 : 
cov(i,r) _ E((i—u,)(r—u,)) a E(ir)-E@E(r) (1) 
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其 中 ，i 和 rr 分别 表示 在 输入 图 像 和 参考 图 像 相同 位 置 的 小 区 域内 的 像素 的 灰 度 值 ， 
这 两 个 区 域 的 形状 称 为 p-kernel, 它 是 一 个 正方 形 。u, Flu, 分别 表示 各 自 小 区 域内 所 
有 像素 的 灰 度 值 的 均值 , o; 和 ea, 分 别 表示 各 自 小 区 域内 所 有 像素 的 灰 度 值 的 标准 差 。 
cov 表示 计算 协 方差 , 已 表 示 计 算数 学 期 望 。 表示 p-kernel 的 像素 数量 。 使 用 式 (2) 
来 计算 皮尔 逊 相关 系数 。 使 用 这 些小 区 域内 的 像素 的 灰 度 值 来 计算 其 中 心 像 素 位 置 
Gor) 的 皮尔 逊 相关 系数 已 ，〈 值 为 -1 到 1) 。 对 于 边界 情况 ， 采 用 截断 处 理 ， 即 忽 
略 边 界 以 外 的 像素 ， 减 少 p-kermel 内 的 像素 数量 n 。 最 后 ， 将 得 到 一 个 与 原始 图 像 尺 
寸 大 小 相同 的 皮尔 逊 相关 系数 矩阵 五 ， 简 称 相关 图 。 一 般 情 况 下 ， 当 皮尔 逊 相关 系 
数 小 于 0 时 ， 表 示 负 相关 ; 当 皮 尔 逊 相关 系数 大 于 0 时 ， 表 示 正 相关 。 皮 和 尔 逊 相关 系 
数 的 绝对 值 越 大 ， 相 关 性 越 强 ， 反 之 越 弱 。 

由 于 想 要 尽 可 能 地 消除 两 幅 图 像 间 相似 的 部 分 ， 并 保留 不 相似 的 部 分 ， 本 文 使 
用 一 个 修改 过 的 sigmoid 函 数 来 修改 皮尔 逊 相关 系数 的 值 。 有 具体 操作 方式 如 下 : 


(3) 
l+e~* 
1 
BR < fee = (4) 


其 中 ， 式 (3) 是 sigmoid 函 数 。 式 (4) 是 修改 后 的 sigmoid 函 数 ， 其 函数 图 象 如 图 2 
红色 实 线 所 示 ， 横 轴 表 示 原 始 的 皮尔 逊 相关 系数 的 值 ， 纵 轴 表 示 修 改 后 的 皮尔 逊 相 
关系 数 的 值 ， 黑 色 虚 线 表 示 没 有 修改 皮尔 逊 相关 系数 的 转换 关系 。 皮 和 尔 逊 相关 系数 
的 值 经 过 修改 后 的 相关 图 如 图 1 中 子 图 (e) Aras, 此 时 相关 图 的 对 比 度 变 得 更 大 了 ， 
有 利于 消除 两 幅 图 像 间 流 量 分 布 相似 的 部 分 ， 并 保留 流量 分 布 不 相似 的 部 分 。 需 要 
注意 的 是 ， 原 始 的 皮尔 逊 相关 系数 可 能 为 负 值 ， 但 其 绝对 值 不 会 很 大 ， 因 为 对 于 两 
幅 已 经 对 齐 、 分 别 减 去 自己 的 背景 和 归 一 化 的 图 像 来 说 基本 上 没有 强 负 相关 。 使 用 
这 种 修改 过 的 sigmoid 函 数 可 以 有 效 地 消除 负 值 的 影响 。 

当 运 动 目标 与 其 他 天 体重 有 登 的 时 候 ， 相 关 性 也 会 出 现 较 强 的 情况 ， 这 并 不 利于 
保留 差异 。 为 了 解决 这 一 问题 ， 引 入 了 差分 流量 系数 。 其 工作 原理 如 下 : 


ee LO (5) 
a min( i, 9r) 
0, dif... =T. 

F, “ee (6) 


RP, i 和 7 分 别 表 示 输 入 图 像 和 参考 图 像 相同 位 置 的 小 区 域内 的 像素 的 灰 度 值 ， 这 
两 个 区 域 的 形状 被 称 为 d-kernel， 它 也 是 一 个 正方 形 。 对 于 边界 处 理 ， 采 用 同 求 皮 尔 
还 相关 系数 一 样 的 原则 。 显 然 ， 可 以 用 式 (5) 来 得 到 两 个 小 区 域内 流量 变化 的 
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(e) 
图 1 (a) 输 入 图 像 ，(b) 参 考 图 像 ，(c) 差 异 图 ; 
(dd) 相关 图 ; (e) 使 用 修改 过 的 sigmoid 函数 修改 后 的 相关 图 ; (ORE R 
Fig.1 (a)input image; (b)reference image; (c)difference image; 
(d)correlation image; (e) correlation image using modified sigmoid function; (f)threshold image 


”图 2 修改 后 的 sigmoid 函数 图 象 

Fig.2 modified sigmoid function 
Ak diff, .例如 , 当 一 个 d-kerel 内 所 有 像素 的 灰 度 值 的 总 和 为 13, 男 一 个 d-kernel 
内 所 有 像素 的 灰 度 值 的 总 和 为 10 时 ， 这 两 个 小 区 域内 流量 变化 的 百分比 为 30%。 在 
式 (6) 中 设置 一 个 阐 值 7〈 立 值 7 的 设置 详 见 3.2 节 ) . Sof, KEETTN, F, 


为 0， 和 否则 ， 瓦 ,为 1。 利 用 这 些小 区 域内 的 像素 的 灰 度 值 来 计算 其 中 心 像素 位 置 
lior) 的 差分 流量 系数 玉 ，。 最 后 ， 将 得 到 一 个 与 原始 图 像 尺寸 大 小 相同 的 差分 流 
m= ASUS A, MMAR. HARA ISA O 所 示 ， 红 色 圆 圈 〈 即 运动 目 
标 ) 内 的 数值 为 0。 当 运动 目标 的 重 双 区域 的 流量 变化 百分比 大 于 等 于 TT 时 ， 阐 值 图 
可 以 用 来 使 运动 目标 在 重 儿 区 域 的 强 相关 失效 , 即 削 弱 相 关 性 , 以 便 突出 运动 天 体 。 

基于 以 上 所 得 的 相关 图 和 阔 值 图 ， 假 设 有 上 剔除 了 坏 像 素 和 饱和 像素 的 输入 图 像 
7 和 参考 图 像 R。 本 文 算法 的 步骤 如 下 : 

步骤 1: I 和 尺 分 别 减 去 自己 的 背景 值 ， 然 后 归 一 化 得 到 了 和 R,。 相 关 的 表达 
式 如 下 : 


I,, =1-bkg(D, (7) 
S, =abs(sum(,,,)), (8) 
I, =o. (9) 

S; 
R,» = R-bkg(R), (10) 
Sp = abs(sum(R,,)), (11) 
R = 各 (12) 

SR 


bkg 表示 求 图 像 的 背景 值 , 本 文 使 用 Bradley 等 人 外 开发 的 Photutils 软 件 求 图像 的 背景 
值 , 其 他 的 方法 也 可 以 尝试 用 来 获取 图 像 的 背景 。 下 标 nb 表 示 去 除 背 景 ，5, 和 Si 表 
示 归 一 化 的 缩放 比例 ，abs 表示 求 绝 对 值 ，sum 表 示 求 图 像 所 有 像素 的 像素 值 的 总 
和 和。 注意， 去 除 坏 像素 和 饱和 像素 的 目的 是 为 了 确保 归 一 化 过 程 的 正确 性 。 

步骤 2: 找 出 1 和 R, 之 间 像 素 灰 度 值 的 差异 D,。 计 算 公式 如 下 : 

Dy =1y—R. (13) 

步骤 3: 消除 两 幅 图 像 间 流量 分 布 相似 的 部 分 ， 并 保留 流量 分 布 不 相似 的 部 分 。 

具体 操作 如 下 : 


P, = pearson(1,, R,), (14) 
F, = diff (1,,R,), (15) 
D,=D,(1-B,f). (16) 
D=D5S,. (17) 


pearson Nit EAA RP. FP Ash (4) feck, dif 表示 计算 两 
HEAVAMT EEA F. A (16) 用 来 消除 两 幅 图 像 之 间 具 有 相似 流量 分 布 的 部 分 ， 并 
保留 具有 不 相似 流量 分 布 的 部 分 ， 且 用 阔 值 图 来 使 运动 目标 在 重合 区 域 的 强 相关 失 
效 。 由 于 在 步骤 1 中 对 万 进行 了 归 一 化 , 需要 对 它 进 行 恢复 , 即 乘 以 系数 9, FBI D COL 
式 (17) ) ， 刀 就 是 本 文 算 法 得 到 的 加 工 后 的 减法 图 像 ， 简 称 esapd 减 法 图 像 。 


2 算法 性 能 

为 了 验证 本 文 算法 的 性 能 ， 我 们 选择 了 三 份 代码 与 本 文 算 法 进行 比较 。 基 于 
Alard 方法 的 代码 有 ISIS2.2(http://www2?.iap.fr/users/alard/package.html) ~ 
diapl(https://users.camk.edu.pl/pych/DIAPL/) 和 HOTPANTS (https://github.com/acbeck 
er/hotpants)， 本 文选 择 HOTPANTS 作为 Alard 算法 的 代表 ， 因 为 它 的 接口 使 用 起 来 
比较 方便 。OIS(https://github.com/quatrope/ois) 中 提供 了 基于 Bramich 方法 的 代码 ， 
而 基于 Zackay 等 人 的 方法 的 代码 是 properimage(https://github.com/quatrope/ProperI 
mage)。 本 文 将 这 三 种 算法 称 为 Alard、Bramich 和 properimage。 需 要 指出 的 是 ,esapd、 
Bramich 和 Alard 都 是 运行 在 2020 MacBook Pro M1 上 。 由 于 安装 环境 的 原因 ， 
properimage 算法 运行 在 台式 计算 机 (Intel(R) 核 心 (TM)i7-10700CPU@2.90GHz) 上 ,其 
运行 时 间 约 为 2020 MacBook Pro M1 的 1.25 倍 。 在 比较 四 种 算法 的 运行 时 间 时 ,将 
对 properimage 算法 的 运行 时 间 乘 以 0.8， 从 而 可 以 相对 准确 地 比较 四 种 算法 的 运行 
时 间 。 接 下 来 通过 试验 比较 展示 本 文 算法 的 性 能 。 
2.1 对 运动 目标 的 检测 

测试 了 2013 年 2 月 5 日 在 云南 天 文 台 丽江 2.4 米 望 远 镜 拍摄 的 Apophis AK. 
由 于 Apophis 的 运动 速度 非常 快 ， 我 们 可 以 用 esapd 图 像 减法 技术 清楚 地 看 到 
Apophis 的 移动 情况 。 本 文 使 用 裁 前 后 的 900x900 的 图 像 进行 试验 .对 于 properimage 
和 HOTPANTS， 使 用 默认 参数 。Bramich 使 用 (21,21) 的 卷 积 核 ，esapd 使 用 p-kernel 
大 小 为 7x7，d-kernel 大 小 为 15x15。 因 为 本 文 的 目的 是 寻找 运动 目标 Apophis, fit 
以 设置 7 为 一 个 较 大 的 值 0.5， 以 确保 尽 可 能 地 消除 较 大 的 残 差 和 去 除 伪 对 象 。 

如 图 3 所 示 ， 应 用 了 esapd 算法 的 减法 图 像 非常 干净 。 在 esapd 减法 图 像 中 ， 白 
点 是 输入 图 像 中 的 运动 目标 Apophis， 黑 点 是 参考 图 像 中 的 运动 目标 Apophis。 在 
esapd 减法 图 像 中 ， 我 们 可 以 很 容易 地 找到 运动 目标 。 然 而 ， 由 算法 properimage、 
Bramich 和 Alard 产生 的 三 幅 减法 图 像 并 不 干净 , 给 寻找 运动 目标 带 来 了 困难 。 通 过 
SAOImageDS9 软件 5 观察 减法 图 像 的 像素 值 发 现 ，esapd 算法 可 以 减 得 非常 干净 ， 
而 其 他 方法 不 能 做 到 这 一 点 。 大 振幅 的 残 差 干扰 了 对 目标 Apophis 的 检测 。Bramich 
和 Alard 引起 的 减法 伪 影 使 得 很 难 确定 一 个 瞬 变 源 是 否 真 实 , 根据 Zackay 等 人 的 研 
究 ， 这 是 由 Alard 与 Lupton 系列 方法 的 不 对 称 性 导致 的 。 在 本 文 算法 中 , 我 们 只 需 
要 使 用 DAOFIND 程序 "来 搜索 运动 目标 Apophis， 而 其 他 三 种 算法 很 难 用 该 程序 
或 其 他 程序 来 搜索 运动 目标 Apophis 。 

找到 的 运动 目标 Apophis 的 位 置 显示 在 esapd 减法 图 像 中 的 红色 圆圈 中 。 值 得 注 
意 的 是 ， 在 参考 图 像 的 右 侧 有 一 个 小 黑 点 (绿色 圆圈 内 ) ， 这 是 由 于 使 用 的 望远镜 
的 CCD 的 缺陷 造成 的 。 事实 上 ， 黑 点 的 值 只 比 图 像 背 景 值 9850 小 200 左右 。 


esapd 
running time: 0.493s 


properimage Bramich Alard 
running time: 5.085s running time: 41.074s $ running time: 1.005s 


3 2013 年 2 月 5 日 拍摄 的 Apophis 图 像 的 图 像 减法 测试 结果 。 从 左 到 右 (顶部 ): 输入 图 像 , 参考 图 像 和 esapd 
减法 图 像 ，( 底 部 )proprimage 减法 图 像 , Bramich 减法 图 像 和 Alard 减法 图 像 

Fig.3 Image subtraction results for test on Apophis images taken on February 5, 2013. Left to right(top): the input 
image,the reference image, and the esapd subtraction image; (bottom)the properimage subtraction image, the 
Bramich subtraction image, and the Alard subtraction image 


通过 人 眼 观察 ， 发 现 只 有 properimage 算法 才能 检测 到 黑 点 ， 因 此 值得 肯定 的 是 ， 
properimage 算法 在 某 些 方面 优 于 本 文 算 法 ， 特 别 是 在 细节 检测 方面 。 但 总 的 来 说 ， 
本 文 算法 具有 良好 的 鲁 棒 性 ， 而 且 本 文 的 算法 运行 速度 最 快 。Bramich 算法 运行 的 
速度 很 慢 。 当 我 们 调整 卷 积 核 的 大 小 为 〈31,31) 时 ， 我 们 需要 运行 191.744 秒 。 

为 了 验证 本 文 的 算法 得 到 的 运动 目标 的 位 置 的 准确 性 ， 我 们 通过 测量 减法 图 像 
中 的 目标 和 输入 图 像 中 的 目标 之 间 的 位 置 偏差 来 进行 试验 。 例如, 测量 图 3 中 esapd 
减法 图 像 的 红色 圆圈 中 的 白 点 位 置 及 其 在 输入 图 像 中 的 对 应 位 置 ， 然 后 计算 它们 的 
位 置 偏差 ， 得 到 本 文 算法 的 位 置 偏差 。 同 理 ， 可 得 到 其 他 算法 的 位 置 偏差 。 本 文 使 
用 2013 年 2 月 4 日 在 云南 天 文 台 丽 江 2.4 米 望 远 镜 拍摄 的 74 幅 Apophis 图 像 进行 
试验 ， 它 们 分 别 来 自 5 个 视 场 。 我 们 按 视 场 将 它们 分 为 5 组 ， 然 后 选取 每 一 组 的 最 
后 一 幅 图 像 作 为 参考 图 像 ， 其 余 的 作为 输入 图 像 。 用 输入 图 像 减 去 参考 图 像 ， 得 到 
减法 图 像 。 使 用 二 维 高 斯 拟 合算 法 02 测 量 运动 目标 Apophis 的 像素 位 置 。 图 4 显示 
了 本 文 得 到 的 结果 , 每 种 算法 都 有 69 个 位 置 偏差 。 令 人 惊讶 的 是 ， 这 四 种 算法 得 出 
的 位 置 偏差 在 均值 、 标 准 差 和 偏 移 量 上 基本 一 致 。 需 要 说 明 的 是 ， 出 现 较 大 位 置 偏 
差 的 原因 是 运动 目标 发 生 了 重 倒 。 本 文 esapd 算法 的 位 置 偏差 在 x HAD y 轴 上 都 有 
最 小 的 标准 差 。 在 x 轴 上 ,本 文 esapd 算法 得 到 的 均值 为 0.0572, 标准 差 为 0.1811 。 
在 y 轴 上 ， 本 文 esapd 算法 得 到 的 均值 为 0.0498， 标 准 差 为 0.1794。 
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图 4 顶部 子 图 : 运动 目标 Apophis 在 减法 图 像 和 输入 图 像 之 间 的 x 轴 位 置 偏差 。 每 种 算法 共有 69 个 位 置 偏差 。 
平均 值 和 标准 差 显示 在 子 图 的 右上 角 。 底 部 子 图 : 与 顶部 子 图 相同 ， 但 是 对 于 y 轴 

Fig.4 Top panel: the object Apophis position difference in x-axis between subtraction image and input image. A total 

of 69 position difference for each algorithm. The mean and standard deviation are shown in the upper right 

comer of the panel. Bottom panel: same as Top panel, but for y-axis 


可 以 看 出 ， 本 文 esapd 算法 得 到 的 目标 具有 相对 稳定 的 位 置 。 
2.2 运行 时 间 对 比 

比较 四 种 算法 的 运行 时 间 。 将 从 ZTF (Zwicky Transient Facility) 下载 的 M35 
星团 图 像 的 大 小 裁剪 为 256x256，512x512，.…..，3072x3072 进行 计算 ， 然 后 得 到 
算法 的 运行 时 间 , 单位 为 秒 。 从 图 5 中 可 以 看 出 , 本 文 esapd 算法 的 运行 时 间 最 快 ， 
最 慢 的 是 Bramich 算法 。 当 图 像 过 大 的 时 候 , Bramich 算法 甚至 会 出 现 无 解 的 情况 。 
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图 5 不 同 尺寸 大 小 的 图 像 在 每 个 算法 的 运行 时 间 (单位 : 秒 ) 


Fig.5 Overview of running time(unit: second) of each algorithm for images of different size 

3 解释 与 讨论 
3.1 背景 ， 归 一 化 和 对 齐 

背景 去 除 和 归 一 化 是 两 个 重要 的 操作 ， 如 果 没 有 这 两 个 操作 ， 本 文 算法 就 无 法 
实现 。 归 一 化 主要 用 于 解决 不 同 曝 光 时 间 图 像 的 问题 。 需 要 注意 的 是 ， 我 们 的 代码 
并 没有 提供 图 像 配 准 算 法 和 背景 估计 算法 。 本文 使 用 Beroiz 等 人 5 开发 的 图 像 对 齐 
算法 astroalign 进行 图 像 的 配 准 工作 和 Bradley 等 人 外 开发 的 Photutils 软件 进行 1-D 
或 2-D 背景 估计 ， 其 他 的 图 像 配 准 算法 和 背景 估计 算法 也 可 以 被 尝试 使 用 。 

需要 注意 的 是 ， 对 齐 精度 对 相关 系数 具有 一 定 影响 。 本 文 使 用 Photutils 生成 一 
幅 大 小 为 14x14, 只 含有 一 颗 星 且 其 中 心 在 图 像 中 心 的 模拟 星象 图 (简称 基准 图 )， 
然后 生成 大 量 类 似 14x14， 但 其 包含 的 星 的 中 心 相 对 于 图 像 中 心 偏 移 了 的 图 像 ， 偏 
移 步 长 (x 和 y 方向 均 可 偏 移 ) 为 0.1 个 像素 (简称 目标 图 ) 。 最 后 利用 这 些 目标 图 
依次 与 基准 图 求 相关 图 ， 并 对 每 个 相关 图 的 所 有 像素 值 求 平均 值 ， 将 此 平均 值 作为 
偏 移 后 的 相关 系数 。 如 图 6 Aras, 例如 ， 当 x 轴 和 y 轴 的 偏 移 量 为 C0, 0) 的 时 候 ， 
相关 系数 最 大 ， 值 为 1; Sx WMM y 轴 的 偏 移 量 为 (2.5，2.5) 的 时 候 ， 相 关系 数 的 
值 约 为 0.85。 实 验 发 现 ， 偏 移 量 越 大 ， 即 对 齐 精度 越 低 ， 相 关系 数 〈 未 使 用 sigmoid 
函数 修改 ) BN; 反之， 相关 系数 越 大 。 知 把 相关 系数 视 为 因 变 量 ， 偏 移 量 视 为 自 
变量 ， 因 变量 与 自 变 量 之 间 类 似 高 斯 分 布 。 这 一 结论 说 明 ， 对 章 精度 越 高 ， 得 到 的 
相关 系数 越 准 确 。 
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图 6 相关 系数 与 偏 移 量 的 关系 


Fig.6 Relationship between correlation coefficient and shift 


3.2 kernel 与 阐 值 7 的 设置 及 算法 的 适用 性 

p-kernel 的 大 小 与 两 幅 图 像 的 半 高 全 宽 (FWHM) LAAR, ZRK, p-kernel 
越 大 ， 反 之 ，p-kernel 越 小 。 例 如 ， 在 同一 晚上 拍摄 的 图 像 中 ，FWHM 相差 不 大 的 
情况 下 ,可 以 设置 p-kernel 为 最 小 值 3x3。d-kernel 与 两 幅 图 像 的 平均 FWHM 有 关 : 
平均 FWHM 越 大 ，d-kernel 越 大 ; 反之 ，d-kernel 越 小 。 

经 实践 发 现 ，p-kernel 以 2-3 倍 FWHM 之 差 ，d-kernel 以 2-3 倍 平均 FWHM 来 
设置 效果 最 佳 , 即 在 突出 运动 目标 的 前 提 下 减 得 干净 。 例如， “4 FWHM 之 差 为 2.36 
pixel 时 ， 可 以 设置 p-kernel 大 小 为 Sx$， 当 平均 FWHM 为 7.08 pixel 时 ， 可 以 设置 
d-kernel 为 15x15。 当 然 ， 此 结论 仅 在 处 理 模 拟 图 像 下 得 出 ， 在 处 理 真实 图 像 时 ， 
还 需 配 合 阔 值 7 来 找到 合适 的 参数 。 

闵 值 7 越 小 ， 突 出 两 幅 图 像 的 差异 的 能 力 越 强 。 然 而 ， 它 并 不 是 越 小 越 好 ， 如 
SET 的 值 太 小 ， 减 法 图 像 就 会 出 现 较 大 振幅 的 残 差 ， 甚 至 是 伪 对 象 ， 因 为 此 时 相关 
性 容易 失效 ， 算 法 退化 为 两 幅 图 像 直 接 相 减 。 通 常 ， 将 7 的 值 设 置 为 0.1。 如 果 两 幅 
图 像 的 差异 较 大 ， 又 想 在 减 得 干净 的 前 提 下 找到 运动 目标 ， 可 将 了 设 大 一 点 ， 比 如 
0.5， 这 里 的 0.5 也 可 以 理解 为 某 一 区 域 的 流量 变化 超过 50% 且 无 论 相关 性 强 弱 如 何 ， 
两 幅 图 像 的 差异 都 会 被 显示 出 来 。 

由 于 本 文 算法 的 特殊 性 ， 当 两 幅 图 像 的 点 扩散 函数 (PSF) 差异 过 大 的 时 候 ， 例 
如 ， 一 个 是 圆 形 的 高 斯 分 布 PSF， 男 一 个 是 倾斜 的 椭圆 高 斯 分 布 PSF 的 时 候 ， 本 文 
算法 的 处 理 效果 并 不 理想 ,但 是 ,如果 两 幅 图 像 都 是 倾斜 角度 一 样 的 椭圆 高 斯 分 布 ， 
本 文 算法 处 理 效果 几乎 不 受 影响 。 本 文 算法 处 理 FWHM 相差 太 大 的 图 像 的 效果 也 
不 是 很 理想 ， 建 议 使 用 FWHM 之 差 在 2 个 像素 之 内 的 图 像 。 

3.3 积分 图 

为 了 加 快 皮 尔 逊 相关 系数 和 差分 流量 系数 的 计算 速度 ,参考 了 Viola 与 Jones 提 出 
的 一 种 利用 积分 图 快速 计算 Haar-like 特 征 的 方法 中。 从 公式 (2) 和 (5) 中， 我 们 
可 以 看 到 有 许多 连 加 运算 ， wmi Drs LPs Lr 和》 ir， 本 文 用 这 种 方法 
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来 加 快 连 加 运算 。 积 分 图 的 使 用 不 仅 提 高 了 皮尔 逊 相关 系数 和 差分 流量 系数 的 计算 
速度 ， 而 且 一 旦 生成 积分 图 ， 当 p-kernel 或 d-kernel 的 大 小 变化 时 ， 运 行 时 间 是 不 会 
因此 改变 的 ， 即 本 文 算法 的 时 间 复 杂 度 与 参数 无 关 。 
3.4 噪声 估计 

减法 图 像 中 的 噪声 有 两 个 来 源 ， 即 输入 图 像 中 的 噪声 和 参考 图 像 中 的 噪声 。 假 
设 两 幅 图 像 的 噪声 符合 泊 松 分 布 c= VN 。 本 文 推导 出 的 减法 图 像 的 噪声 估计 公式 


如 | : 
Rp 


注意 ， 该 公式 仅 适 用 于 已 ,= 0 或 相关 强度 很 弱 (皮尔 逊 相关 系数 的 值 小 于 0.2) 的 情 
况 。 
4 总 结 与 展望 

本 文 提 出 了 一 种 可 供 选择 的 图 像 减法 算法 esapd, 期 望 将 其 应 用 于 搜索 运动 目标 。 
该 算法 旨 在 利用 相关 性 消除 两 幅 图 像 间 具 有 相似 流量 分 布 的 部 分 ， 保 留 具有 不 相似 
流量 分 布 的 部 分 。 试 验 表明 ,该 算法 可 以 在 很 短 的 时 间 内 找 出 两 幅 图 像 之 间 的 差异 ， 
并 搜索 出 运动 目标 。 同 时 ， 该 算法 存在 一 些 不 足 ， 如 果 两 幅 图 像 的 PSF 分 布 差异 过 
大 ， 本 文 算法 处 理 的 效果 有 待 提高 ， 因 此 ， 后 续 工作 将 结合 使 用 PSF 来 进一步 研究 
天 文 图 像 减 法 技术 。 
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Abstract: In the field of astronomical image processing, image subtraction techniques 
are widely used. Due to the varying seeing conditions, different exposure times, and 
different filters, the two images cannot do image subtraction directly. A novel image 
subtraction algorithm is proposed in the paper. Its essence is based on correlation to 
eliminate the parts with the similar flux distribution and preserve the parts with the 
different flux distribution between two images in sense of statistics. The algorithm can be 
fast executive, is numerically stable, and is locally independent. Based on our algorithm, 
we have developed a set of image subtraction code with Python as the interface and C as 
the implementation. We use the algorithm, together with three other similar algorithms, 
to perform with the astronomical images . Experiments show that our algorithm can find 
the difference between two images and detect moving objects in a very short time while 
having good robustness and position measurement stability. 
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